import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 定义原始函数
def original_function(x):
    return np.sin(x*5) * x-np.random.normal(0, 0.1, x.shape)
#np.exp(-x*np.sin(x))+
# 定义圆谐函数拟合函数
def circular_harmonics(theta, *coeffs):
    """
    圆谐函数拟合函数
    :param theta: 角度变量
    :param coeffs: 拟合系数，coeffs[0] 是直流分量，coeffs[1::2] 是 cos(kθ) 的系数，coeffs[2::2] 是 sin(kθ) 的系数
    :return: 拟合值
    """
    result = np.ones_like(theta)*coeffs[0]  # 直流分量
    for k in range(1, len(coeffs) // 2 + 1):
        result += coeffs[k] * np.cos(k * theta) + coeffs[k + len(coeffs) // 2] * np.sin(k * theta)
    return result

def decompose(x,y,num_harmonics=25):
    initial_guess = [0] * (2 * num_harmonics + 1)  # 初始猜测系数
    coeffs, _ = curve_fit(circular_harmonics, x, y, p0=initial_guess)
    return coeffs

import csv
import math
def test():
    x=np.linspace(0,2*np.pi,600+1)
    x=x[:-1]
    y=np.ones_like(x)
    z=np.zeros_like(x)
    y_n=np.cos(x)
    y_t=-np.sin(x)

    dx = 49.994499200000000
    dy = 0.52441693000000000
    dr=math.sqrt(dx*dx+dy*dy)
    t=math.asin(dy/dr)

    k1=y_n*np.cos(x)+y_t*np.cos(x+np.pi/2)
    k2=y_n*np.sin(x)+y_t*np.sin(x+np.pi/2)
        # 保存 (y_n, z) 到 CSV 文件
    with open('yn_z.csv', 'w', newline='') as file:
        writer = csv.writer(file)

        for i in range(len(y_n)):
            writer.writerow([y_n[i], z[i]])

    # 保存 (z, y_t) 到 CSV 文件
    with open('z_yt.csv', 'w', newline='') as file:
        writer = csv.writer(file)
    
        for i in range(len(z)):
            writer.writerow([z[i], y_t[i]])

    # 保存 (y_n, y_t) 到 CSV 文件
    with open('yn_yt.csv', 'w', newline='') as file:
        writer = csv.writer(file)
        
        for i in range(len(y_n)):
            writer.writerow([y_n[i], y_t[i]])

    
    #载荷分解
    coeffs=decompose(x,y_n,8)
    y_fit=circular_harmonics(x,*coeffs)
    print(coeffs)

    coeffs=decompose(x,y_t,8)
    y_fit=circular_harmonics(x,*coeffs)
    print(coeffs)


    n=4
    for i in range(1,n+1):
        print(i)
        k = (i - 1) * 1.0 / n 
        print(np.pi * 2.0 * k)

def test_cmp():
    #i=1...600
    angle=[]
    for i in range(1,11,1):
        k=(i-1)/10
        angle.append(np.pi*2.0*k)
    print(angle)
    x=np.linspace(0,2*np.pi,10+1)
    print(x)

  # 运行主函数
if __name__ == "__main__":
    test()
